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Quantum Monte Carlo methods have proven to be valuable in the study of strongly 
correlated quantum systems, particularly nuclear physics and cold atomic gases. Historically, 
such ab initio simulations have been used to study properties of light nuclei, including spectra 
and form factors, low-energy scattering, and high-momentum properties including inclusive 
scattering and one- and two-body momentum distributions. More recently they have been 
used to study the properties of homogeneous and inhomogeneous neutron matter and cold 
atomic gases. There are close analogies between these seemingly diverse systems, including 
the equation of state, superfluid pairing, and linear response to external probes. In this 
paper, we compare and contrast results found in nuclear and cold atom physics. We show 
updated lattice results for the energy of the homogeneous unitary Fermi gas and comparisons 
with neutron matter, as well as for the dependence of the cold atom energy on the mass ratio 
between paired particles, which yields insights on the structure of the ground state. We also 
provide new lattice and continuum results for the harmonically trapped unitary gas, again 
comparing neutron matter and cold atoms. 

§1. Introduction 

Quantum Monte Carlo methodic 'PI have been proven very valuable in study- 
ing a host of strongly correlated quantum systems, including solid and liquid He- 
humjlf'0'[!f'§ electronic systems|3l§ light nuclei)®© and more recently neutron 
matteJUJ'lllS'lll'Ill and cold atomic Fermi gases 113 'Hll'lllUM Many of these stud- 
ies are formulated in the continuum to be able to describe the short-range repulsion 
between, for example, helium atoms, nucleons, or electrons. For low-density neu- 
tron matter and cold Fermi gases, where the dominant interaction is a relatively 
short-range attraction, lattice methods have also proven valuableI13.|llf'[lll'lllf In 
this paper we describe applications of Quantum Monte Carlo methods to the equa- 
tion of state, superfluid pairing gap, and related properties of neutron matter and 
cold atomic gases. 

We concentrate on the zero temperature properties of cold atoms and neutron 
matter. Cold Fermi atoms have a deceptively simple, essentially zero-range interac- 
tion, H = —h/2m\/f -\- Y^ - ■ Vo5{Ti — Tj), the strength of which (Vb) can be tuned 
to produce a very rich set of physics described by a relatively small set of universal 
parameters. As the strength of the (attractive) interaction Vq increases, one goes 
from a weak BCS pairing regime to the BEC regime of strongly bound pairs. Many 
simulations and experimentaHSf are performed near and at the unitary limit, which 
is where the two-body system produces a nearly zero-energy bound state. 
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The s-wave interaction between two neutrons is also very attractive, nearly pro- 
ducing a bound state. For very dilute neutron matter, the neutron matter and cold 
atom equations-of-state should be very similar as a function of the product of the 
fermi momentum kp = (Svr^p)^/^ times the magnitude a of the scattering length. 
The s-wave interaction between neutrons also has a significant effective range which 
causes the equation of state of neutron matter and cold atoms to diverge at relatively 
modest densities. The effective range of the interaction also affects the pairing gap, 
etc. It may be possible to use narrow Feshbach resonances to more directly mimic 
neutron matter and study the dependence of the equation of state on the effective 
rangeHH* There are, of course, additional p— wave interactions between neutrons, 
though these are relatively modest at low densities. 

§2. Monte Carlo Methods 

Monte Carlo methods have proven quite effective in dealing with strongly cor- 
related quantum systems. They are most efficient when they incorporate as much 
knowledge of the physical system to be studied as possible. The zero-temperature 
Diffusion Monte Carlo (DMC) and Auxiliary Field Monte Carlo (AFMC) methods 
we employ use Monte Carlo to propagate a trial wave function to the true ground 
state of quantum systems: 

N 

|l^o> = ^M-Ht] \'^t) = n ^M-H{t/N)] \^t), (1) 

i=l 

where the imaginary time propagation is split into small imaginary time steps 5t = 
t/N. For small 5t the propagator can be evaluated accurately in terms of the two- 
body propagator: 

(R' I ex.p[-H6T] I R) 

9iji^ijT'''ij) 

where gij and g^j are the interacting and free two-particle propagators respectively, 
determined by the eigenstates of the interacting and free two-particle Hamiltonian: 
hij = + v{rij) and h^^ = — ^Vf^-. The factor exp[— /i?(5r] is the free single particle 
propagator, a simple gaussian, and R is a 3A^-dimensional vector containing the 
coordinate-space positions of all the particles. 

Diffusion Monte Carlo and Auxiliary Field Monte Carlo operate in different 
spaces: DMC performs the simulations in coordinate space, Monte Carlo methods 
are used to sample the spatial integrals. AFMC calculations are carried out on a 
lattice, and can be viewed as the evolution of single-particle orbitals in imaginary 
time. In these lattice calculations the single-particle orbitals can be transformed 
between coordinate and momentum space through fast fourier transforms. In AFMC 
the Monte Carlo is incorporated by evaluating the potential matrix elements in terms 
of fluctuating auxiliary fieldsHi 



(r'l exp[-/iij()>] |r). 
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In either case, it is possible to use a very accurate short-time two-body propaga- 
tor. In the continuum DMC simulations and in the limit of zero-range interactions 
the propagator at unitarity takes the simple formlH^ 



5.,(r:,, r,,) = 4.(4-, r,,) + ^ , ^ eM-m/ {nHr){n,^ + r^')/4], (3) 

where the correction to the free particle propagator arises from the pair propagating 
to the same point and then diffusion plus multiple scattering terms. In actual simu- 
lations typically an analytic potential is used and a short-time approximation used 
for the propagator: 

5i,(r^ r,,) = exp[-F(r^)5r/2] 4.(4, r,,) exp[-y(r,,)5r/2]. (4) 

This expression is accurate to order ((5r)^. 

For lattice calculations, the simplest propagator used is for an on-site attractive 
interaction and either a Hubbard-like hopping Hamiltonian or a k'^/{2m) kinetic 
term in an expression analogous to Eq. [4j 

gij{r'ij,rij) = exp[-r5r/2] exp[-V6T] exp[-r5T/2]. (5) 

The potential is evaluated by Monte Carlo sampling of auxiliary fields and the ki- 
netic energy exactly through the use of fast fourier transforms IHll This simplified 
interaction yields a finite range of the order of the lattice spacing. One can remove 
this residual effective range by altering either the kinetic term T, adding higher-order 
momentum terms in the kinetic energy, or by introducing additional auxiliary fields 
to give the correct low-energy two-body spectrumJUf'HII' Lattice methods cannot be 
fully galilean invariant as there is a lattice cutoff at the high- momentum scale, and 
simple implementations of the improved actions can have similar effects to a finite 
effective range for pairs with non-zero momentumllU These corrections typically 
vanish in the limit of large lattices. 

The AFMC lattice simulations have the great advantage that they do not suffer 
from a sign problem for purely attractive two-body interactions. The up and down 
spin evolution can be factorized in AFMC and for unpolarized systems the spin 
up and spin down determinants are real and equal and the product is positive. 
Therefore, the Monte Carlo results should be exact within statistical errors if we 
can perform simulations that are sufficiently dilute, for sufficiently large number of 
particles, and at low-enough temperatures (large r). 

The DMC continuum simulations, in contrast, suffer from a sign problem as the 
spin up and spin down determinants are independent and the product can be positive 
or negative. The overlap of the propagated wave function with the trial function has 
a statistical error that grows with respect to the propagation time and the number 
of particles. The fixed-node algorithm, though, provides an accurate upper bound 
to the energy by requiring the paths in the evolution not to cross planes where the 
trial wave function is zero. This can be proved to provide an upper bound to the 
true ground state energy, and hence the trial wave function can be optimized as the 
one that produces the lowest energy. The advantage of the DMC simulations is that 
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the fixed-node algorithm is equahy applicable to polarized or unpolarized systems, 
to systems with unequal masses, etc., where the AFMC method also suffers from a 
sign problem. Similar constrained-path algorithms exist for AFMC methodsjl^J'llll' 
but do not typically provide upper bounds. 

In both DMC and AFMC simulations it is very important to use good trial 
functions in strongly correlated ground state calculations. In DMC the fixed-node 
results are quite accurate if a good trial function of BCS type is employed, as we 
shall discuss later. In AFMC the results should be correct for any trial function that 
has a finite overlap with the true ground state. However the statistical errors are 
dramatically reduced if one uses a BCS trial function for the trial state, as described 



in Ref . 23 ) . It is also extremely valuable to use a branching random walk algorithm 
to limit statistical errors in both the DMC simulations, as traditionally done, and 
in the AFMC simulationsHSl'IlD In zero-temperature calculations we are trying to 
reach the eigenstate of the transfer matrix exp[— ifr], the branching random walk 
algorithm is a Markov chain algorithm, that is, only the most recent history of the 
path is required for performing the next step. Consequently one can iterate to very 
low temperatures/large imaginary times compared to other algorithms. 

For quantities other than the energy, we often evaluate matrix elements of the 
form: 

{^T\Oexp[-HTm 

= {^r\eM-Hrm) ' 
where iZ'j is an initial state used to start the simulation and \I^t is a trial state 
incorporating as much knowledge as possible of the ground state. It is often possible 
to take 'I^i = though this is not required. In the limit of large imaginary time r, 
the 0{t) can be used to determine the ground-state properties of the system. The 
energy is the simplest: in that case the Hamiltonian commutes with the propagator 
exp[— ffr] and one can calculate the ground state expectation value. For other 
properties one must insert the propagation symmetrically between the initial and 
final state or create a new Hamiltonian H' = H + eO, in either case evaluating the 
expectation value of the operator: 

which in the limit of large r and t' gives the true ground state expectation value. 



§3. Equation of State 



The equations of state for cold atoms and for neutron matter have been exten- 
sively studied theoretically, and the cold atom system has been extensively studied 
experimentally as well. The cold atom system is very simple: for a zero-range 
interaction the equation of state is a function of only the product of the Fermi 
momentum and the scattering length (kpa). More specifically, the energy can 
be written as a function of kpa times the energy of a non- interacting Fermi Gas 
{Efg = {2./b){h^/2m)kl) at the same density. Here the Fermi momentum is de- 
fined through the density of the corresponding non-interacting gas: kp = (Svr^/))^/^. 
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At unitarity (infinite scattering length) the ratio of the energies of interactin g a nd 

1321 



noninteracting Fermi gases E/Epc is typicaUy called the Bertsch parameter ^ 
3.1. Unitarity 




Fig. 1. AFMC lattice calculations of the unitary Fermi Gas f parameter, updated from Ref. |23[ ). 
Symbols are for different kinetic terms as a function of particle number and lattice size. The 
lattice spacing is denoted as a. Simulations have been performed with lattices, for different 
values of lattice length L in each direction; open symbols are for even L=16,20,24; closed are 
for odd L (see text). All extrapolations are consistent with ^ = 0.372(5). 



A history of results for the Bertsch parameter is given in Ref. 28). The first 
DMC calculation used up to 40 particles and a modified Poeschl- Teller potential 
with kpTe ~ 0.3, where is the effective range of the interaction, and yielded a 
fixed- node energy of ^ = 0.44(1) 03 Subsequent DMC calculations used improved 
trial functions, larger particle numbers, and better extrapolations to kpTe — t- to 
yield ^ = 0.40(1)1^ The best present DMC result is from the calculations of Ref. 



34), while an updated extrapolation to rg — )■ gives ^ = 0.390(1)2^ for an upper 
bound. This calculation also carefully compared results at finite particle number to 
a superfluid Local Density Approximation (LDA) to extrapolate to large N. It was 
found that calculations for = 38 or larger are very close to the thermodynamic 
limit. 

There is also a substantial history of lattice simulations, both for the ground- 
stateJlll''[llf'llZf'lll}'[MI and at finite temperatureMJ'Ulf The earliest ground-state cal- 
culations estimated = 0.25(3), for systems up to 22 particles on lattices up to 



6^. The recent calculations of Ref. 23l) use branching random walks and a BCS 



trial function and importance sampling for systems of 66 particles on lattices up to 
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27^ and obtain ^ = 0.372(5) for several different actions. Updated results for these 
calculations are shown in Figure [T] 

In the figure, the upper curves use a k'^ dispersion relation tuned to unitarity. 
This k'^ dispersion has a finite positive effective range of 0.337 a, where a is the lattice 
spacing. The middle set of curves adopt a, k"^ + k'^ dispersion that is tuned to zero 
effective range, and the lower curves use a Hubbard dispersion relation, which has a 
negative effective range of —0.306 a. The k"^ + A;^ results show a set of simulations 
with even L as open symbols, while simulations at odd L are shown as filled symbols. 
The two sets of results are slightly displaced; similar displacements have been found 
with limited statistics for the other dispersions. All extrapolate to the same value 
of within statistical errors; we return to the dependence on effective range below. 
A new lattice calculation in Ref. 28) reports a higher value of ^, above the upper 
bound found in DMC calculations. 

There have also been a large number of experimental determinations of ^: the 
original measurementdl^'ED'EH' have found qualitative agreement with the DMC 
calculations listed above. More precise recent experiments have found ^ = 0.39(2jlll 
and ^ = 0.41(ljlil with a smaller value of ^ = 0.375(5) found most recentlylllf This 
experimental value is quite precise and overlaps our lattice results. 

3.2. Equation of State: Cold Atoms and Neutron Matter 

Of course the full equation of state {E/Efg) as a function of kpa is required to 
compare with neutron matter, which has a fixed, large effective range and must be 
studied by varying the density. The most recent DMC results for the full equation 
of state are presented in Fig. [2j and compared to the lattice results and the most 
recent experimental result. These results are quite smooth as a function of kpa and 
extrapolate correctly in both the BCS and BEC regimes. 

Because the cold-atom interaction is short-ranged, the derivative of the energy 
with respect to kpa is given completely by short-range physics, as originally written 
down by Tan in a series of papers®'©'!!!! The derivative of the energy per particle 
with respect to kpa is given, using the Hellman-Feynman theorem, by: 

dE N f . , dV(r) 

The pair distribution g^i{r) — )• goes like A? jr^ at unitarity for small r, with 
g^\,{r) — )• 1/2 at large r. The change in energy with respect to is 

dE _ h^27,pA^ .^_«„2„2.2 



da ^ m 



C = Stt'p'A^ (9) 



where C is Tan's contact parameter. Near unitarity the EOS is conventionally 
parametrized as 

^ e-T^ + ..., (10) 



EpG kpa 

with C = {5tt/2)C /kp. We return to the contact parameter in the discussion of 
short-range physics below. 
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Fig. 2. Equation of state of cold atoms versus l/ikpa). Blue circles are DMC calculations, the red 
square and green diamond are lattice and experimental values at unitarity l/{kFa) = 0. The 
insert shows the corrections from finite effective range near unitarity (see text). 



In Fig. [sjthese cold atom results are compared to the QMC for neutron matter jlll 
and to the analytic expression available at small kpa. At low densities, the neutron 
matter and cold atom results agree, they also agree with a simple extrapolation of 
the analytic results near kpa = 0. At higher densities, the cold atom and neutron 
matter equations of state start to diverge somewhat as the effective range becomes 
important. The dependence of the equation of state on effective range can be made 
explicit, as we discuss below. This dependence gives a quantitative picture of the 
difference between neutron matter and cold atoms that could perhaps be tested in 
cold atom experiments with narrow resonances. We will return to the finite-range 
corrections below. 

3.3. Equation of State: Unequal Masses 

Cold atom experiments can also be performed with species of different mass, 
providing important information about the structure of the ground state of the 
unitary Fermi Gas. For species of different mass and mi, if we normalize the ^ 

parameter by the reduced mass Epc = ~^jf~' BCS theory would give a value of ^ 
independent of the mass ratio r = m^/mi. The difference in Hamiltonians for equal 
(r = 1 , m-|- = mi = m) masses and unequal masses is 

^ 2mi ^ 2mi ^ 2m 
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Fig. 3. Comparison of the equation of state of cold atoms and neutron matter at low density. 



Neutron matter calculations are from Ref. 14 1. Differences at low density are primarily due to 
the effective range of the neutron-neutron interaction. The solid line is a fit to the cold atom 
results, the dashed line includes an estimate of effective range effects (see text). 

#pairs 2 /-„ iN2 

= V (11) 

where in the last line the particles have arbitrarily been divided into N/2 spin up - 
spin down pairs. 

Figure |4] shows the DMC calculations of ^ for different mass ratios. Initial 
calculations for different mass ratios were reported in Ref. |49) ). From Eq. Ill we can 



see that the energy change can be evaluated in perturbation theory near r=l. 

A{E/N) = {Q\AH\ 0) = {l/2){PfJ{Am))\rJ^^^^, (12) 



where the 1/2 comes from the number of pairs (N/2), Pij is the total momentum of 
a pair, and the expectation value is to be taken in the ground state of the equal mass 
unitary gas. Note that these calculations were performed for small but finite value 
of the effective range, yielding a slightly larger value of ^ than at zero effective range. 
This rather asymmetric way of writing the energy difference is valuable because it 
tells us something about the character of the state. For a BCS-like state with all 
pairs at P = the energy difference is zero in first-order perturbation theory. Of 
course the free Fermi Gas can also be written in this manner. The difference is finite 
for the case when the ground state wave function does not have a spin down particle 
at — p for every spin up particle at momentum p. 



9 



0.4 

0.395^^ 
0.39 
"J" 0.385 
0.38 

0.375 
0.37, 



0.1 0.2 0.3 0.4 

(r-l)'/(r+lf 



0.5 



0.6 



Fig. 4. Dependence of the unitary Fermi Gas equation of state on mass ratio for fixed reduced mass 
/i = mfmi/{m^ + mi), plotted versus (r — Vf /{r + 1)^, where r is the mass ratio. 

Fig. |4] shows the DMC calculations as points with error bars, and a quadratic 
fit to the data. The linear coefficient in this fit is very small, consistent with zero 
within statistical errors. Thus to a very good approximation the ground state of the 
unitary gas can be written as a state of pairs with zero momentum. To confirm this 
result it would be important to have experimental measurements of the energy for 
several different mass ratios. 

3.4. Equation of State: effective range 

As is apparent from Fig. [T| the equation of state for finite effective range 
varies linearly with kpVf, at small effective range. In Ref. 23^, the equation of state 
at unitarity for different effective ranges was found to be: 



iikpre) = m + Skpre 



(13) 



where ^(0) characterizes the ground-state energy at zero effective range, and S is 
the slope parameter giving the linear dependence on kpre- The slope parameter was 
extracted from both DMC and AFMC calculations^) and found to be 5 = 0.12(3). 
The results for different effective ranges are shown in Figure [5} More recent DMC 
results for a variety of interactions have recently appeared,^'' ' they find S = 0.127(4) 
using a variety of interactions. These calculations further demonstrate that S is a 



universal parameter, as originally conjectured in the original version of Ref. 29) 



In this manuscript, the authors also claim that lattice results will in general have 
a dependence on the total momentum P of a pair. For the unitary gas, however, 
the expectation value of (-P^) is approximately zero, as shown above. Therefore the 
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lattice and continuum results are both in agreement with Eq. [T3| 
0.42 




Fig. 5. Dependence of the unitary Fermi Gas equation of state on Fermi momentum vs. effective 
range {kpre)- Sliaded bands are fits to the lattice results, and dashed lines give DMC results. 

In DMC calculations the slope parameter S is not too sensitive to kpa near 
unitarity. Fig. [2] shows, in the inset, the slope parameter S evaluated from DMC 
calculations near unitarity. It is positive and approximately 0.1 near unitarity, but 
changes significantly in the BCS and BEC regimes. The difference between the cold 
atom EOS and neutron matter at sufficiently small densities should be approximately 
^neutrons — (.atoms ~ Skpr^, or approximately 0.05 at —kpa = 5 since the neutron- 
neutron effective range is expected to be approximately 2.7 fm. Fig. [3] shows a fit 
to the cold atom results at zero effective range as a solid line. The dashed line adds 
an effective range correction with S = 0.1. This should be the dominant correction 
at kp < 0.25 fm~^ , near kp ~ 0.5 fm~^ one would have kprg ~ 1 and higher order 
corrections in s— and p— wave interactions could be important. 



§4. Pairing Gap 



Both low-density neutron matter and cold atoms are strongly paired Fermi sys- 
tems, they exhibit some of the largest pairing gaps of any systems known when 
measured in terms of the Fermi energy. We define the pairing gap at T=0 as the 
difference between the energy of an odd particle system and the average of the two 
nearby even particle systems in periodic boundary conditions: 

A = EiN + 1) - {E{N) + E{N + 2))/2, (14) 

with the universal parameter 5 defined as the pairing gap divided by the Fermi energy 
Ep = h?kp/2m. For simulations of a large enough number of particles this should 
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correspond to the traditional definition of tlie pairing gap. 

Tlie fact that the pairing gap is so large, a sizable fraction of the Fermi energy, 
makes it possible to use QMC methods to accurately calculate the gap by separately 
calculating the energies of the even and odd particle systems. In addition, the fact 
that the energy per particle shows no significant shell effects for reasonably small 
systems {N > 30) makes it much easier to approach the continuum limit. Though 
there is an upper bound principle for the even and odd systems, there is no specific 
bound on the pairing gap. 

The original calculations of the pairing gap in cold atoms at unitarity found 
A/EpQ « 0.9 oi 5 = 0.55(5)113* Subsequent improvements to the wave functiorHD 
found a slightly reduced value for the gap, 5 = 0.50(5). These results can be com- 
pared to an extraction of the pairing gap from the measured density distributions in 
partially spin-polarized trapped cold atoms^'^ ' and measurements of the RF response 
in such systems jm who find 6 = 0.45(5) and 5 = 0.44(3), respectively. 

The pairing gap in neutron matter has historically been the subject of a great 
deal of interest and theoretical activityMf'llll' QMC calculations of the pairing gap 
were performed in 13) and|14[). These calculations used the s-wave and s- & p-wave 
components o f the AV18 interaction, respectively. A summary of the results are 
shown in Fig. pll2f 




Fig. 6. Pairing gap in cold atoms and neutron matter. BCS mean-field results are shown as solid 
lines, DMC results are shown as symbols. The pairing gaps are divided with the relevant one- 
body quantity, namely the Fermi energy Ep (analogously to the ground-state energy in Fig. [TT] 
being divided with Efg ~ SEf/5. 

In the figure, BCS results are given by solid lines. In the weak-coupling limit, 
the pairing gap is expected to be reduced from the BCS value by (l/4e)^/'^ 0.45 
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from the Gorkov polarization correctional' It is difficult for QMC calculations to 
calculate the pairing gap at coupling weaker than kpa = —1 because of the delicate 
cancellations. At this coupling, though, we find a suppression in the gap roughly 
compatible with the Gorkov suppression. At stronger coupling the suppression di- 
minishes smoothly, and the gap reaches a value of 0.50(05) at unitarity. In the BEG 
regime the pairing gap approaches half the binding energy of the pair as reproduced 
by the BGS equation. 

The calculated gaps in neutron matter are considerably smaller than in cold 
atoms, but still reach a maximum of nearly O.SEp at kpa = —5. The effective 
range in the neutron-neutron interaction reduces the gap significantly, as shown in 
the comparison of BGS results and in the Monte Garlo calculations. These pairing 
gaps are considerably larger than found in many diagrammatic approaches but 



in agreement with the lattice results of Ref. 24 ) 



§5. Short-range physics 



One can also investigate the short-range (high-momentum) physics in cold atom 
experiments. Because of the simple short-range interaction, this short range physics 
is directly related to the equation of state discussed above. The pair distribution 
function and off-diagonal density matrix at short distances and the momentum dis- 
tribution at high momenta are all governed by the contact parameter. 

The pair distribution function is shown in Fig. [7} the contact governs the huge 
spike near r = 0. In the inset the pair distribution function is multiplied by (kpr)'^ 
to show the behavior near r = 0, the dip at very short distances is due to the finite 
range interaction used in the simulations. The behavior in the region of the vertical 
dashed line and slightly beyond is governed by universal physics. In the figure the 
red points are variational Monte Garlo (VMG) results obtained from the trial wave 
function \^t), the green are the mixed estimates of the form obtained from Eq. [6] 
and blue are the full extrapolated DMG results, as is apparent in the figure the VMG 
calculation is correctly capturing the basic physics. 

The momentum distribution scaled by /c^ is plotted in Fig. [s} The momentum 
distribution at large k is proportional to the contact. The horizontal line in the 
figure is the value that would be expected from calculations of the equation of state: 
C = 0.901(2). Extractions of the contact from all these observables are consistent 
with this value, though of course some are noisier than others. It will be interesting to 
see what information on the short-range physics of neutron matter can be obtained 
from theory and experiments with narrow resonances with a significant effective 
range. 

Initial experiments on the spin and density response of cold atoms have also 
been performedHll These and future results will be very interesting as they can tell 
us about the propagation of particles and spin in the unitary gas. The response 
functions can be written as: 



S,iq,u;) = (0| J]exp[-iq • r,]|/ )(/| J]exp[^q • r,-]|0) 6{u; - {Ef - Eq)) 

f i j 
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Fig. 7. Pair distribution function g^i{r) for cold atoms at unitarity. Inset shows the behavior at 
short distances scaled by r^, the magnitude of the contact determines the value of this quantity 
(see text). 

S^{q, u;) =Y, (0| exp[-iq • r^] cr, |/ ) • (/| ^ exp[iq • r,] aj |0) S{u; - {Ef - Eq)) 

f i 3 

(15) 

These response functions have been calculated at high momenta in terms of the 
operator production and related high- momentum expansions The 
experiments show a two-peak structure in the density response, one at w = / (2m) 
associated with the breaking of a pair and one at w = /{Am) associated with the 
propagation of a pair. The spin response requires breaking of a pair. The initial 
experiments are at rather high momentum transfer, many times the Fermi momen- 
tum, and hence probe the short-range physics. It will be interesting to see how these 
response functions evolve at lower momenta. One can also compare calculations and 
experiments on the sum rules associated with the contact parameter by integrating 
the response over u. 

§6. Inhomogeneous Matter 

Finally, we turn to the properties of inhomogeneous matter. This is the regime 
with perhaps the closest connection between nuclear physics and cold atom physics. 
In the inner crust of a neutron star the neutrons form a gas between the neutron- 
rich nuclei that exist in a lattice structure. This neutron matter is very low density 
and is inhomogeneous, and for many of the transport properties we would like to 
understand the behavior of this gas. The properties of inhomogeneous neutron mat- 
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Fig. 8. Momentum distribution scaled by A;* for cold atoms at unitarity. 



ter, particularly the gradient terms, are very difficult to determine from the binding 
energies of atomic nuclei. The isovector gradient term is one of the least constrained 
parameters in nuclear density functionals, and ab-initio calculations can provide 
valuable guidance. 

6.1. Inhomogeneous Matter: Bulk Properties 

We have calculated the properties of ffiiite systems of neutrons bound in har- 



monic and Woods-Saxon wellsJHf Original calculations of these dropa^ were limited 
to N=8 neutrons because of the spin-dependence of the nuclear interaction. These 
more recent calculations use Auxiliary Field Diffusion Monte Carl(B@'[i6j methods 
and Green's function Monte Carlo methods; the former have been used to treat up 
to N=50 neutrons. The AV18 NN interaction plus the UIX three-nucleon interaction 
have been used for these calculations. As these wells produce fairly modest densities 
for 50 particles or less, the three-neutron interaction is not very important. 

The energies in Fig.jojhave been scaled by 1/(^A^^/^), the expected behavior in 
the Thomas-Fermi approximation for an EOS of the form E/Epc = The upper 
horizontal dashed lines are for ^ = 1 (free fermions), and the lower for ^ = 0.5. 
Neutron matter over a considerable range of densities is roughly consistent with 
^ = 0.5, though it is less attractive at low and high densities. The traps have 
harmonic frequencies of a; = 10 MeV (upper panel) and 5 MeV ( lower panel). 
Results for several typical older-generation density functionals are shown as open 
circles, and calculations using SLY4 as solid lines. 

The Skyrme interactions typically give significantly lower energies than the mi- 
croscopic calculations, particularly for the 10 MeV well. The greater difference for 
the 10 MeV well suggests that isovector gradient terms in the density functional 
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Fig. 9. Scaled energies for neutrons bound in a harmonic wellH) The upper and lower straight 
dashed lines are the Thomas-Fermi mean-field results for free fermions and for a scaled EOS 
E/N = ^Efg with ^ — 0.5. The open symbols are calculations with previous generation density 
functionals, and the filled symbols are GFMC and AFDMC calculations. The jagged line shows 
the results for the SLyi density functional. The upper panel is for a harmonic trap frequency 
of 10 MeV and the lower for 5 MeV. 



should be more repulsive. The curves marked "SLY4-adj" in the figure are obtained 
by adjusting the isovector gradient, the isovector pairing, and the isovector spin- 
orbit terms in the interaction. At the closed shells {N = 8, 20, 40) only the change 
in the isovector gradient term is important. A reasonable fit can be obtained to the 
closed-shell energies in the 5 and 10 MeV harmonic wells and the Woods-Saxon wells 
with a single adjustment to the isovector gradient term. This adjustment also better 
reproduces the rms radii and mass distributions of the ab-initio calculations!^ 

It is interesting to compare the neutron drop results to those obtained for cold 
atoms. In Fig. 10 we plot the scaled energies of cold atom systems obtained in DMC 
and AFMC calculations as well as the results of previous calculations£Z<''i5£<''iHZi' Two 
previous DMC calculations used fairly simple trial wave functionsjl3'[l§ the first 
using an orbital basis for the trial function and the second using a simple 1/r BCS 
pairing function as the trial wave function. Both resulted in energies far above what 
would be expected in local-density (Thomas Fermi) theory with ^ ~ 0.4. New lattice 
result^ yield somewhat lower energies, but significant shell structure. They are 
also above the energies expected from the measured and calculated values of ^ in the 
continuum, indicating either unnaturally large gradient corrections or other effects. 

Our new DMC and AFMC calculations produce energies considerably lower than 
previous results. The DMC energies are very smooth as a function of N as compared 
to neutron drop results, indicating a lack of shell closures for the unitary Fermi gas 
with zero effective range. These DMC calculations use a more sophisticated trial 
function incorporating both a single-particle shell model basis and substantial short- 
range pairing into the trial wave function. The resulting energy is somewhat higher 
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than the expectations from the bulk for N up to 50 particles. 

It has been showrIZ3 that the cold atoms trapped in harmonic wells have a 
breathing mode associated with the scale invariance of the Hamiltonian of exactly 
2 huj, independent of particle number N, giving further evidence that cold atoms at 
unitarity have a very weak shell structure, if any. Initial AFMC lattice calculations 
give energies for 30 fermions much closer to the expected bulk limit, indicating a 
smooth and rapid transition from few-particle systems to the bulk. Further DMC 
and AFMC calculations are being pursued. 
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Fig. 10. Comparison of different calculations of the harmonically confined unitary Fermi 
GasHlUllUlll'lllUlHl The vertical axis is the square of the energy divided by Thomas-Fermi 
energy Etf = iu{3N)^^^ /4, the horizontal axis is the number of particles. For a very large 
system (E/Etf)^ should approach ^. 



In Fig. [TT] we compare cold atom results to neutron drops in the same traps 
as a function of particle number A^. We plot the square of the energies because, 
in the local density approximation, the square of the energy of the confined system 
is proportional to the energy of uniform matter determined by The bulk limit 
as obtained from the lattice calculations shown in Fig. [T] is shown as an arrow at 
the lower right of the figure. The energies of neutron drops are considerably higher 
than those of cold atoms, this is at least partially the result of the effective range 
in neutron matter. The gradient terms are likely also important, however. These 
can be more precisely constrained by performing a local density calculation using a 
realistic equation of state for neutron matter. 

6.2. Inhomogeneous Matter: pairing 

We have also performed calculations of the pairing gap in neutron drops, shown 



in Fig. 12 In atomic nuclei pairing is, at least predominantly, a bulk effect, in that 



the coherence length appears to be comparable to the size of the nucleus. This is to 
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further corrections. 



be expected in a regime where the effective range is comparable to the interparticle 
spacing, in such cases it should be possible to construct a mean-field theory that 
produces a qualitatively correct picture. 

For cold atoms the conclusions will be quite different; for large enough systems 
the unpaired atom will necessarily sit outside the center of the drop. The gap 
is simply too high in the high-density central region for the unpaired particle to 
penetrate. It will be very instructive to compare theories and experiments as a 
function of scattering length, effective range, and particle number. 



§7. Conclusions 



Cold atom experiments and theory provide many valuable insights into our un- 
derstanding of strongly correlated fermions, and in particular have a close relation- 
ship with low-density neutron matter. The equations of state for low-density neutron 
matter and cold atoms are by now well understood and very similar as a function 
of kpa, and the difference is understood in terms of a correction proportional to the 
fermi momentum times the effective range. 

Pairing gaps in cold atoms demonstrate a smooth transition from the BEC to 
BCS regime, and indicate that a sizable pairing gap is to be expected. Further 
experimental and theoretical studies as a function of kpTe at and near unitarity would 
be very valuable in providing explicit confirmation and a more precise understanding. 

A better understanding of the linear response of cold Fermi atoms could lead to 
new insights into the dynamic response of neutron matter. It would be particularly 
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Fig. 12. Calculated pairing gaps in harmonically confined neutron drops, at harmonic frequencies 
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closure. Gaps are calculated from the mass difference formula with a phase factor of —1^+^ to 
take into the account the fact that the unpaired particles have a higher energy associated with 
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valuable to map out both the density and spin response of cold atoms as a function 
of momentum transfer. Analogies to topics such as neutrino propagation in dense 
matter are clear, though there is not a direct correspondence as there is for the 
equation of state and pairing gap. 

Inhomogeneous matter is also quite intriguing, including small systems of trapped 
fermions, fermions in optical lattices, and the transition from three to two-dimensions. 
Inhomogeneous cold atom systems also have close analogies in nuclear physics, in- 
cluding the physics of nuclei and the neutron star crust. Studies of narrow resonances 
with finite effective range could help us understand the evolution of pairing from a 
local to a bulk phenomenon. The rapidly expanding scope of cold atom experiments 
and Quantum Monte Carlo will undoubtedly reveal intriguing new physics and close 
correlations with nuclear physics. 
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